#---Load Packages----
library(lavaan)
library(xtable)
setwd("~/Replication Code/")

#---Loading Data----
anes <- read.dta("./data/anes_subset.dta")
anes16 <- read.dta("./data/anes16_subset.dta")

#---Mode Model----
invar_dat <- subset(anes16, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "f2f"))

FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "f2f"

# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc
# reverse wording
thard ~~ specfavr
'
# dless seems best anchor given loading
fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit


fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("pdisc ~ 1"),
                       group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar_part, standardized = T)
round(fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
scalar_fit_part <- fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit_part <- rbind(comb_fit, scalar_fit_part)
comb_fit_part

# Combining results for output
FIT_INVAR <- c("chisq","df", "cfi", "srmr", "rmsea")
invar_fit <- data.frame(rbind(fitMeasures(fit_config)[FIT_INVAR],
                              fitMeasures(fit_metric)[FIT_INVAR],
                              fitMeasures(fit_scalar)[FIT_INVAR],
                              fitMeasures(fit_scalar_part)[FIT_INVAR]))
invar_fit$chi2_change <- NA
invar_fit$chi2_change_p <- NA
invar_fit$cfi_change <- NA
invar_fit$srmr_change <- NA
invar_fit$rmsea_change <- NA

invar_fit$chi2_change[2] <- invar_fit$chisq[2] - invar_fit$chisq[1]
invar_fit$chi2_change[3] <- invar_fit$chisq[3] - invar_fit$chisq[2]
invar_fit$chi2_change[4] <- invar_fit$chisq[4] - invar_fit$chisq[2]
invar_fit$chi2_change_p[2] <- max(pchisq(invar_fit$chisq[2] - invar_fit$chisq[1], 
                                         invar_fit$df[2] - invar_fit$df[1], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[3] <- max(pchisq(invar_fit$chisq[3] - invar_fit$chisq[2], 
                                         invar_fit$df[3] - invar_fit$df[2], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[4] <- max(pchisq(invar_fit$chisq[4] - invar_fit$chisq[2], 
                                         invar_fit$df[4] - invar_fit$df[2], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[2] <- paste0("$<$", invar_fit$chi2_change_p[2])
invar_fit$chi2_change_p[3] <- paste0("$<$", invar_fit$chi2_change_p[3])
invar_fit$chi2_change_p[4] <- paste0("$<$", invar_fit$chi2_change_p[4])
invar_fit$cfi_change[2] <- invar_fit$cfi[2] - invar_fit$cfi[1]
invar_fit$cfi_change[3] <- invar_fit$cfi[3] - invar_fit$cfi[2]
invar_fit$cfi_change[4] <- invar_fit$cfi[4] - invar_fit$cfi[2]
invar_fit$srmr_change[2] <- invar_fit$srmr[2] - invar_fit$srmr[1]
invar_fit$srmr_change[3] <- invar_fit$srmr[3] - invar_fit$srmr[2]
invar_fit$srmr_change[4] <- invar_fit$srmr[4] - invar_fit$srmr[2]
invar_fit$rmsea_change[2] <- invar_fit$rmsea[2] - invar_fit$rmsea[1]
invar_fit$rmsea_change[3] <- invar_fit$rmsea[3] - invar_fit$rmsea[2]
invar_fit$rmsea_change[4] <- invar_fit$rmsea[4] - invar_fit$rmsea[2]


# removing DF
invar_table <- invar_fit[,-2]
# formatting to generate table
colnames(invar_table) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", "$\\Delta\\chi^2$", "$\\Delta\\chi^2$ p-value", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")
rownames(invar_table) <- c("Configural", "Metric", "Scalar", "Scalar (Partial)")

tab <- xtable(invar_table, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment Over Time",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_year_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "s", "fg", "fg", "fg"))
# NOTE <- list()
# NOTE$pos <- list()
# NOTE$pos[[1]] <- c(nrow(invar_table))
# NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{past discrimination} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item intercepts. Data from the 1992 and 2016 ANES. One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      hline.after = c(-1, 0))

#---Temporal Model----
invar_dat <- subset(anes, 
                    year == 1992, 
                    select = c("specfavr", "pdisc", "thard", "dless", "year"))
invar_dat16 <- subset(anes16, f2f == 1, 
                      select = c("specfavr", "pdisc", "thard", "dless"))
invar_dat16$year <- 2016
invar_dat <- rbind(invar_dat, invar_dat16)

FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")


GROUP <- "year"
mod <- '
rr =~ pdisc + dless + thard + specfavr
# reverse wording
thard ~~ specfavr # seems appropriate
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]


fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus",
                  
                  group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

fit_scalar_part <- cfa(mod, invar_dat, mimic = "Mplus", 
                       group.partial = c("specfavr ~ 1", "thard ~ 1"), # for 1992-94; 1992-2016
                       group = GROUP, group.equal = c("loadings", "intercepts"))
summary(fit_scalar_part, standardized = T)
round(fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")], 3)
scalar_fit_part <- fitMeasures(fit_scalar_part)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit_part <- rbind(comb_fit, scalar_fit_part)
comb_fit_part


# Combining results for output
FIT_INVAR <- c("chisq","df", "cfi", "srmr", "rmsea")
invar_fit <- data.frame(rbind(fitMeasures(fit_config)[FIT_INVAR],
                              fitMeasures(fit_metric)[FIT_INVAR],
                              fitMeasures(fit_scalar)[FIT_INVAR],
                              fitMeasures(fit_scalar_part)[FIT_INVAR]))
invar_fit$chi2_change <- NA
invar_fit$chi2_change_p <- NA
invar_fit$cfi_change <- NA
invar_fit$srmr_change <- NA
invar_fit$rmsea_change <- NA

invar_fit$chi2_change[2] <- invar_fit$chisq[2] - invar_fit$chisq[1]
invar_fit$chi2_change[3] <- invar_fit$chisq[3] - invar_fit$chisq[2]
invar_fit$chi2_change[4] <- invar_fit$chisq[4] - invar_fit$chisq[2]
invar_fit$chi2_change_p[2] <- max(pchisq(invar_fit$chisq[2] - invar_fit$chisq[1], 
                                         invar_fit$df[2] - invar_fit$df[1], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[3] <- max(pchisq(invar_fit$chisq[3] - invar_fit$chisq[2], 
                                         invar_fit$df[3] - invar_fit$df[2], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[4] <- max(pchisq(invar_fit$chisq[4] - invar_fit$chisq[2], 
                                         invar_fit$df[4] - invar_fit$df[2], 
                                         lower.tail = F), .001)
invar_fit$chi2_change_p[2] <- paste0("$<$", invar_fit$chi2_change_p[2])
invar_fit$chi2_change_p[3] <- paste0("$<$", invar_fit$chi2_change_p[3])
invar_fit$chi2_change_p[4] <- paste0("$<$", invar_fit$chi2_change_p[4])
invar_fit$cfi_change[2] <- invar_fit$cfi[2] - invar_fit$cfi[1]
invar_fit$cfi_change[3] <- invar_fit$cfi[3] - invar_fit$cfi[2]
invar_fit$cfi_change[4] <- invar_fit$cfi[4] - invar_fit$cfi[2]
invar_fit$srmr_change[2] <- invar_fit$srmr[2] - invar_fit$srmr[1]
invar_fit$srmr_change[3] <- invar_fit$srmr[3] - invar_fit$srmr[2]
invar_fit$srmr_change[4] <- invar_fit$srmr[4] - invar_fit$srmr[2]
invar_fit$rmsea_change[2] <- invar_fit$rmsea[2] - invar_fit$rmsea[1]
invar_fit$rmsea_change[3] <- invar_fit$rmsea[3] - invar_fit$rmsea[2]
invar_fit$rmsea_change[4] <- invar_fit$rmsea[4] - invar_fit$rmsea[2]


# removing DF
invar_table <- invar_fit[,-2]
# formatting to generate table
colnames(invar_table) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", "$\\Delta\\chi^2$", "$\\Delta\\chi^2$ p-value", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")
rownames(invar_table) <- c("Configural", "Metric", "Scalar", "Scalar (Partial)")

tab <- xtable(invar_table, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment Over Time",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_year_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "s", "fg", "fg", "fg"))
# NOTE <- list()
# NOTE$pos <- list()
# NOTE$pos[[1]] <- c(nrow(invar_table))
# NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{past discrimination} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains item intercepts. Data from the 1992 and 2016 ANES. One residual correlation estimated between \\textit{try hard} and \\textit{special favors}.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      comment = F,
      hline.after = c(-1, 0))
